Image processing as state reconstruction in optics 
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The image reconstruction of partially coherent light is interpreted as the quan- 
tum state reconstruction. The efficient method based on maximum-likelihood 
estimation is proposed to acquire information from registered intensity mea- 
surements affected by noise. The connection with totally incoherent image 
restoration is pointed out. The feasibility of the method is demonstrated nu- 
merically. Spatial and correlation details significantly smaller than the diffrac- 
tion limit are revealed in the reconstructed pattern. 
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1. Introduction 

Light conveys considerable part of information on the 
surrounding world. Information coded into and trans- 
mitted by means of light plays the key role in contem- 
porary information technology. That is why any deeper 
understanding of fundamental limitations imposed by the 
theory represents a challenging problem. The origin of 
the diffraction limit restricting the spatial resolution is 
comprehended since the era of wave optics. This phe- 
nomenon, manifested for example by a fuzzy diffraction 
spot as a consequence of the finite aperture, yields a se- 
rious limitation in image processing. But there are still 
other physical restrictions which must be taken into ac- 
count. Direct observations are not able to determine the 
phase of optical fields due to the fast oscillations and 
due to the effect of the time integration of intensity de- 
tectors. The phase must be therefore retrieved adopting 
sophisticated techniques. This is known phase problem, 
solution of which is sensitive to noise and requires various 
regularization treatments. 

There are several ways to surpass the mentioned 
limitations imposed by realistic aspects of opti- 
cal observations. Sophisticated algorithms of data 
processing and image reconstruction have been de- 
vised, such as analytic continuation of signal and 
diagonalization of optical dgyiceiiSi&iliiSiLSai, dif- 
ferent regularization and smoothing techniques of 
direct deconvolution and methods of projections 
onto convex sets (POCS) 10 ! 11 ! 12 ! 13 ! 14 ! 15 ! 16 ! 17 ! 18 ! 19 ! 20 , 
solving of the transport-of-intensity equation for 
phase reconstructioriSiiSSiSiiSi 2 ^, utilizing of canonical 
transforms 2 ^, and statistical methods based on minimum 
least-squares distanc o 27 i 28 , maximum entrop}! 2 ^ 3 ^* 3 ^, 
maximum Cramer-Rao bound 3 ^ 3 -, and maximum like- 
lihood principlo 34 i 35 i 36 i 37 i 38 i 39 i 40 i 41 . The tomographical 
synthesis of different intensity observations of an object 
can considerably improve its reconstruction and allow 
the phase retrieval 42 i 43 i 44 i 45 i 46 i 47 i 48 i 49 i 50 i 51 i 52 i 53 . The 
resolution can also be enhanced by eliminating of out-of- 
focus light by means of the apodisation technique 5 ^ or 
by utilizing of confocal arrangement and interaction be- 



tween light and matter, like in multi-photon fluorescence 
microscopy and STED technique 55-56 . 

In this article the problem of image processing will be 
addressed from the viewpoint of statistical reconstruc- 
tion techniques. The series of tomographical-like inten- 
sity measurements will be used up for the reconstruction 
of a state of partially coherent light. In the particular 
case of totally incoherent light the proposed approach 
will be identified with the Richardson algorithm^! of the 
image reconstruction. 

There is a tight connection between fundamental prin- 
ciples of wave optics and quantum mechanics. Descrip- 
tion of the scalar wave in optics is equivalent to descrip- 
tion of the de Broglie wave of a mass particle in the frame- 
work of quantum mechanics. The pure quantum state in 
position representation coincides with the complex scalar 
wave, and similarly, the mixed quantum state in this rep- 
resentation corresponds to the correlation function of the 
second order. This connection will be systematically ex- 
ploited and the problem of quantum state reconstruction 
will be considered in analogy with optical counterpart of 
image processing. 

For the sake of simplicity all the considered problems 
will be treated as two dimensional problems. The first 
dimension corresponds to the evolution parameter — time 
t in dynamical problems or longitudinal z-coordinate in 
the case of image processing. The second dimension cor- 
responds to the observed quantity. This could be position 
in the former case of dynamical problems or transverse 
x-coordinate in the later case of image processing. Fur- 
ther generalization to higher dimension can be obtained 
by a straightforward manner. 

2. Wave theory 

In this section the analogy between scalar wave optics 
and quantum mechanics will be highlighted. As will be 
shown, abstract quantum formulation is advantageous for 
the purpose of signal reconstruction. 

In quantum domain, pure quantum state |^>) from the 
Hilbert space represents the complete knowledge about 
the position and momentum of a particle, of course obey- 
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ing the Heisenberg uncertainty principle. Any random- 
ness in the ensemble of identically prepared particles is 
described by the incoherent mixture of pure states — a 
density operator, 



P 



= ^2^k\fk)('Pk 



(i) 



Probabilities are nonnegative and adds to unity, and 
the mixed state Q satisfies the following relations, 



p, Tr[p] = 1, Mp1V> > 0, V|V), 



(2) 



where + means the Hermitian conjugation. Denoting 
formally the position by a projector complex 
amplitude ip(x) = (x\ip) describes the coherent quasi- 
monochromatic scalar field. Indeed, it characterizes the 
amplitude as well as the phase of the propagating wave. 
In the general case of partially coherent field, the second 
order correlation function^ 

T(x,x') = Tr[p\x')(x\]=J2Mx\fk)(fk\x') = 

k 

= ^2^k(fik(x)(p* k (x / ) = (<ip*(x')ip(x)) ens , (3) 

k 

describes the statistical properties of the field. The 
brackets ( ) ons denote the averaging over complex am- 
plitudes of all modes k. The analogy between the den- 
sity matrix Jy and the mutual intensity J3J is expressed 
clearly by the relations analogous to (J2J, 

r*(x',x) = F(x,x'), Jdxl(x) = l, T(x,x) > 0. (4) 

Notice that the function I(x) — T(x, x) means the optical 
intensity of field. The analogy between quantum and 
wave descriptions can be emphasized in phase space by 
means of the Wigner quasi-distribution£L£&, 



W{x,p) = i Jdx'e- npx 'T{x + x',x-x'). 



(5) 



This (x,p) distribution is real bounded function, which 
is however, not positively defined in general. 

Let us proceed further to consider the state transfor- 
mation. Assuming linearity and causality the equation 
for evolution of pure state formally reads 



|V>)out = T\ip) ia . 
Here T is linear operator satisfying the equation 
8 



dz 



LT, 



(6) 



(7) 



where z is evolution parameter and the generator L of 
evolution is considered to be z-independent. The evolu- 
tion equation JJJ covers both the Schrodinger equation 
in quantum mechanics and paraxial Helmholtz equation 



in Fresnel approximation of scalar wave optics. The uni- 
tary evolution is governed by the Hamiltonian operator, 
L = —j-H, and the evolution of the mixed state is de- 
scribed by transformation 



Pout 



fp in f+, f 



exp 



h 



-Hz 



(8) 



Evolution |JSJ of the state p in the Schrodinger picture 
can be equivalently formulated in the Heisenberg picture. 
This formulation follows the laws of ray optics. Indeed, 
the relationship between the canonical observables of po- 
sition and momentum reads 



•Eout 
Pout 



T 



T. 



(9) 



Particularly, for the evolution generated by quadratic 
Hamiltonian in canonical observables the transformation 
(JSJ is linear, 



•^out 
Pout 



A BUSi, 

C D J I Pi n 



(10) 



where Det[T] = AD - BC = 1. This generic ABCD 
transformation covers useful cases of wave transforma- 
tion, for example free evolution, displacement, rotation, 
phase shift, squeezing, and chirp. Linear transformation 
of the (x, p) operators induces the evolution of the Wigner 
function by linear transformation of its variables, 



W(x,p) = W(Dx -Bp- s, -Cx + Ap 



(11) 



Roughly speaking, all these transformations rotate and 
rescale the input state. As the consequence, the observ- 
able Ax + Bp can be measured offering an important tool 
for all the tomographical methods. 

In the classical limit there is a tight connection be- 
tween evolution of position and momentum operators in 
the Heisenberg picture (|10|l and geometrical paraxial op- 
tics represented by the identity between operators (x,p) 
and its c- values (x,p). The state vector (x,p) is used 
to specify position and angle of the optical ray. Simi- 
lar description may be adopted for particle in classical 
mechanics. Evolution operator T is given by the ABCD 
matrix T completed by the transverse shift s and rotation 
r in analogy with the relation (|10|l . Geometrical optics 
as well as classical mechanics do not involve interference, 
what simplifies the input-output relations considerably. 
This is why the geometrical optics is so suitable for "di- 
rect" observations. Indeed, if the positions x\, x-i for the 
two values z%, z% are measured, the state vector (xo,Po) 
for z = can be completely reconstructed as 



1 



Z2 ~ Zi 



XiZ 2 
X2 



X 2 Zi 
Xl 



(12) 



In the case of losses the evolution turns out to be non- 
unitary. Let us imagine the absorbing screen with 2a 
aperture. The incident state can be decomposed in the 
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base of eigenstates |£) of position x on screen. Since only 
a part of eigenstates spectrum for eigenvalues £ € [—a, a] 
is transmitted, the non-unitary transformation can be 
described by projection operator 



T 



(13) 



corresponding to finite aperture. 

Let us conclude the overview by explicit formulation of 
the state transformation in position representation. The 
generic evolution © of the pure state (coherent wave) 
takes a well-known form of the superposition integral, 



4>out(x) = {x\ 



dx (x|T|x )(x |^)iE 



= j dx h(x,xo)tpin(xo). (14) 

The kernel of the integral transformation l|14|) , 

h(x,x Q ) = (x\T\x ), (15) 

is called propagator in the quantum theory and the re- 
sponse function or point-spread function (PSF) in the 
scalar wave theory^. Loosely speaking, it relates the 
point source in the object (input) plane, z = 0, with 
its image in the image (output) plane with coordinate 
z. This mapping is fuzzy in realistic image processing 
due to the effect of diffraction caused by non-unitary 
evolution. In the case of unlimited aperture it may be- 
come sharp corresponding to the case of unitary evolu- 
tion, T + = T _1 . Analogously, the evolution (JSJ of the 
mixed state (mutual intensity) in the position represen- 
tation reads 

T out (x,x') = J J "dqdq 1 h(x,q)h*(x' ,q')T in (q,q'). (16) 
3. Detection 

According to standard formulation of quantum mechan- 
ics the measurement is represented by an observable, a 
Hermitian operator A. Eigenvalues of this operator cor- 
respond to possible results of elementary measurements. 
Eigenstates determine the possible states after the mea- 
surement and they are complete and orthogonal, 



Ala) = 



^2\a)(a\ = l, (a\a') = S a 



(17) 



These properties are reflected in the probability p a — 
Tr[p\a)(a\] predicted by quantum theory guaranteeing 
the normalization of probabilities ^2 a p a = 1 (complete- 
ness), and mutual exclusivity of the results a (orthog- 
onality). This description may be further generalized 
in terms of positive-operator valued measure (POVM) 
yielding the decomposition of identity operator^SiSi, 



IL, > 0, 



i. 



(18) 



It predicts the probability for registering the output b 
analogously to the case of orthogonal projectors, pb = 
Tr[pIIh]. The notion of POVM plays the crucial role 
in description a generic quantum measurement in state 
estimation and discrimination. 

Registration of the image intensity I{x) of partially 
coherent wave in the transverse position x corresponds 
to the measurement of position operator x in the output 
plane, 



I(x) = Y out (x,x) =p(x) = Tr[p out \x)(x[ 



(19) 



Realistic detector always possesses the finite spatial re- 
solving power. Denoting its pixels by the indices i, the 
simplest representation of detector POVM is given by the 
operators 



Oi = 



A, 



dx \x)(x\, 



(20) 



where the integration is done along the surface of the i-th 
pixel. 

Consider now the generic observation scheme. The in- 
put state p represented by its mutual intensity T(x, x') 
in wave description is transformed by optical device T = 
T(A, B, . . .) with the response function h(x, xq; A, B , . . .). 
Resulting output state p ut is observed by the detector 
(|2*U|l placed in the output plane. The detector counts 
the elementary clicks in every i-th pixel. The numbers 
Ni of registered clicks represented by relative frequen- 
cies fi = Ni/N, N = ^2iNi, sample the probabilities pi 
(intensities /j), 



p % = Tr /SoutA = Tr pTI 



IL = T+0,T. 



(21) 



Notice however, that this scheme is rather classical and 
it does not take into account statistics of detection pro- 
cess in accordance with classical image processing, when 
intensity is considered as measurable quantity. Provided 
that quantum nature of detection will be considered, the 
relation (|20|l should be modified taking into account reg- 
istration of photons instead. 

4. Direct reconstruction 

The reconstruction of the signal in wave theory is rather 
involved and extensive field with many applications. Let 
us review briefly this topic. Assuming an unknown signal 
propagating through optical refractive and diffractive el- 
ements, the output field may be detected. Provided that 
properties of the optical device are known, and detection 
is ideal, the input signal may be predicted from the out- 
put one. This is the classical inverse problem of wave 
optics. 

Standard methods use the isoplanatic approximation 
involving the relation l|14|) as convolution, 



dx h(x - X ) 1pin(x ) = h * Vii 



(22) 
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Inversion is given by the Fourier deconvolution 



Vwt + Af 



(23) 



Here ipi n , Vwt, and h are Fourier transformations of ipi n , 
-0 out , and h, respectively, and A/" represents the spectrum 
of additive noise Af. Typical point-spread function h 
has the form of Sine or BeSinc function and h corre- 
sponds to step-function. Hence, the spatial frequencies 
of the signal are transmitted only up to certain upper 
cut-off 59 . This is why the deconvolution is very sensi- 
tive to noise Af and diverges at spatial frequencies for 
which the transfer function h turns to be zero. In par- 
ticular, for frequencies above the cut-off, the transfer 
function h vanishes and i1-?-SI) diverges due to the broad 
noise spectrum Af. Some regularization procedures are 
necessary in all these cases^ii*i2iiiiiiiSii&. The spe- 
cial attention has been devoted to the more accurate de- 
scription of the optical device. Detailed analysis needs 
special choice of eigenfunctions related to finite aper- 
ture instead of spatial frequencies^* 4 * 5 -. Systematic the- 
ory of this remarkable basis, so called prolate spheroidal 
wave functions, was given by Frieden 6 . Several fur- 
ther super-resolution techniques as apodisation 54 or an- 
alytic continuation^* 8 - have been suggested. The inverse 
source problem may be further generalized, taking into 
account other realistic aspects of detection. For exam- 
ple, optical intensity is detected by real photo-detectors 
instead of complex amplitude and phase is subject of 
reconstruction 1 ! 2 ! 9 ! 16 ! 18 ! 19 ! 20 ! 21 ! 22 ! 23 ! 24 ! 25 ! 26 . 

Standard image processing deals with the observa- 
tion in the image plane revealing the sharpest image. 
However, the observations in defocused planes are also 
worthwhile^* 4 ^* 5 ^* 6 ^. They correspond to the observation 
of Ax + Bp operators in the language of quantum the- 
ory HI U| l , bringing other piece of information about signal 
and affording better employment of measured data. Such 
tomographical technique was suggested by Bertrand, 
Bertrand 4 ^ and Vogel, Riskeni 3 - and experimentally ver- 
ified by group from University of Oregon^. Up to now 
several other tomographical schemes for observing of var- 
ious faces of the system have been proposed. This is usu- 
ally achieved by adjusting of some parameters in the set- 
up. Particular configuration in dependence on parameter 
provides the desired group of transformations. Both the 
classical X-ray tomography (CT) used in medicine 63 ' 64,65 
and the homodyne tomography 4 ^ 4 ^* 5 ^ 5 ^ 5 ^ use the group 
of rotations. Phase-space tomography and chronocyclic 
tomography are connected with the symplectic group but 
they are convertible to the classical tomographji 4 ^* 4 ^* 5 ^. 
General non-homogeneous symplectic tomography was 
introduced by Mancini, Man'ko, and Tombesi-S. 

All the above mentioned inverse problems are relat- 
ing measured data fi with theoretical probabilities pi by 
means of the equality 



where multi-index i passes over all configurations of op- 
tical device and over all pixels of detector. However, the 
solution of linear equations (|24[1 represents an ill-posed 
problem of the same kind as image reconstruction us- 
ing deconvolution. This procedure is very sensitive to 
noise, which is an inevitably involved in any measurement 
scheme. Ill posed problem implies, that reconstructed 
"state" need not represent any physically possible ob- 
ject. In the language of quantum theory this means that 
(ip\p\tfj) < may hold for some states. Alternatively, 
the optical intensity I(x) may drop below zero at some 
position coordinates in wave-theory language. Linear al- 
gorithm is not capable to guarantee such necessary condi- 
tions as positive definiteness of density matrix or mutual 
intensity. 

Statistical approaches suggest a remedy to this prob- 
lem releasing the strict condition (|24|) . For example, the 
equality between /, and pi could be replaced by require- 
ment of their minimal least-squares distance 



El*- 1 * 



(25) 



an obvious choice in engineering practice. But 
the other metrics are also eligible. Least-squares 
methodSL 2 ^ 6 ., Richardson method 3 ^, maximum- 
likelihood (ML) principle and expectation-maximization 
(EM) algorithm 35 ! 36 ! 37 ! 38 ! 39 ! 40 ! 41 ! 67 ! 68 , principle of maxi- 
mum Cramer-Rao bound 3 ^ 3 -, maximum entropy (ME) 
method 2 ^* 3 ^ 3 ^* 6 ^ 7 ^* 7 ^, and intrinsic correlation function 
(ICF) model represent several examples of various 
statistical signal-processing schemes. In the following 
section the arguments in favor of maximum likelihood 
estimation will be formulated with the help of quantum 
treatment. 

5. Reconstruction as generalized quantum mea- 
surement 

Let us assume the generic scheme for the quantum mea- 
surement (|2U p -l]2i p described above. The conditional 
probability of detecting JVj clicks in the i-th pixel, if the 
state p occurs in the input plane, has the form of multi- 
nomial distribution 



£(P)«]> 



Nfi 



(26) 



where fi are the relative frequencies of registered clicks, 
Nfi — Ni. The input state p is the subject of estimation 
procedure. The likelihood functional (|26[) then gives the 
answer to the question "How is it likely that the given 
data fi were registered provided that the system was in 
the given quantum state jo?" For some states the detec- 
tion of given data is more likely than for others. Using 
the relation 121[l the log-likclihood function reads 



Tr pU t 



fii 



(24) 



In C(p) = ^2 fi ln Pi = E f* ln Tr 



(27) 
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Maximum likelihood principle selects such a state p os t for 
which the likelihood reaches its maximum, 



Pest = arg 



max In C(p) 
p 



The formal necessary condition 



5 In C(p) 



6p 



(28) 



(29) 



Pest 



may be rewritten to the form of extremal operator equa- 
^ or> 72 i 73 4 74 i 75 i 76^ or alternatively, extremization can be 
done by means of numerical up-hill simplex method' 7 . 
Any density matrix may be parameterized in diagonal 
form (JIJ using independent (orthogonal) basis states 
\(fk) and the variation l|29l) may be done along these 
rays. Likelihood function depends on the density matrix 
through probabilities pi. This yields the system of cou- 



pled equations 
k. Using the relation 



8 ln£(p) 
<5 (Vk\ 



for any allowed component 



5 ln£(p) 



^ Pi 



(30) 



and the normalization Tr[p] = 1, the extremal 
equationSiZiiZ^iSiZli for the density operator p reads 



Rp = p. 



Here 



R 



~ Pi 



(31) 



(32) 



and probabilities pt are state dependent 1)21(1. Operator 
equation (|31|l determines the most likely solution p es t> for 
which R(p cst ) = 1 holds on the Hilbert space of the state 
Pes J 2]75 - No prior knowledge about the estimated state 
is needed. Results of the measurement itself are sufficient 
for analysis. 

Let us develop the optical counterpart of this recon- 
struction problem. In the spatial domain the extremal 
equation (|31|l has the form of integral equation for mu- 
tual intensity r(x,a/), 



dx"R,(q,x')r(x',q')=r(q,<f), 



(33) 



where the resolution of identity 1 = J dx \x)(x\ has been 
used. Kernel 



TZ(q,x')=^2-V t (q,x') 
~ Pi 



(34) 



and functions 



V i (q,x') = / dxh*(x,q)h(x,x') (35) 



correspond to the operator R and to the POVM operators 
n^, respectively. The probabilities l|21|) of elementary 
detection in individual pixels then read 



Pi = JJdq dq'T(q,q')V l (q', 



(36) 



The equation (|33|l relates measured data /,-, properties 
of optical device, and reconstructed signal T(x,x'). De- 
pendence on the optical apparatus is expressed via point- 
spread function h(x,x') only. However, this mutual re- 
lation is inseparable, since the relation is nonlinear. In 
comparison to standard treatments in scalar optics, no 
assumptions about statistical nature of the signal have 
been done. This seems to be reasonable, since the coher- 
ence properties of the light field may change during the 
propagation (van Cittert-Zernike effect) . The proposed 
formulation anticipates only the knowledge of the opti- 
cal apparatus and the measured data without any prior 
assumptions about the unknown signal. 

Special cases of the generic formulation deserve atten- 
tion. Let us assume an iterative solution of the equa- 
tion l|33|) taking the maximally-ignorant initial guess rep- 
resented by the totally mixed uniform state, p^ ^ — 
jyl, where 1/D ensures the trace normalization in D- 
dimensional Hilbert space. After evaluating the kernel 
we are able to write down the first iteration of es- 
timated state, = _R(°)p<°) = J2i /A /Tr[fy. In the 
spatial domain this state has the form of partially coher- 
ent superposition of response functions (|15fl weighted by 
measured data /i, ^ /, = 1, 



r (1 W) = E/< 



J A dxh*(x,q)h(x,q') 
JdtJ.dx \h(x,0\ 2 



(37) 



It is clear that the coherence properties of estimated sig- 
nal are changed during repeated iterations of extremal 
equation (|33() . Besides the proposed iterative solution the 
well-known EM algorithm^ 7 " completed by unitary sterjSi 
can also be utilized. This guarantees the convergence and 
keeps all fundamental properties l@J of partially coherent 
signal r(a;, x 1 ). 

As the second special case the totally incoherent light 
can be assumed, 



r(x,x') = I(x)S(x-x'), 



(38) 



where 6(x) is Dirac distribution. The extremal equation 
then reduces to 



fdq'K(q,q')I(q')=I(q), 
whereas the probabilities l|36|) read 



dqVi{q, q)H0)- 



(39) 



(40) 



The relations l|39|) . 1)340. and (|40|l provide the extremal 
equations for unknown optical intensity I(x), 

fi 



E 



A, 



SdqVi(q,q)I(q) 



dq'V i (x,q , )I(q')=I(x). (41) 
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This relationship may be utilized for iterative procedure 
as was proposed by Richardson in 1972 for incoherent im- 
age reconstruction. Notice, however, that in the original 
derivation. 34 the Bayes rule was adopted. The treatment 
devised here makes it possible to extend the solution to 
the cases of partially or totally coherent signals. 

6. Numerical example 

In this section we demonstrate the feasibility and advan- 
tages of the presented method by means of the carefully 
selected example. The partially coherent testing object 
is chosen below the resolution limit of the simple optical 
device with finite aperture. Therefore, the observed im- 
ages do not bear any resemblance with the true object. 
In spite of this obvious limitations the reconstructed ob- 
ject reveals the original structure. This improvement is 
achieved by adjusting of the detector position in trans- 
verse and longitudinal directions. The background noise 
is added to simulated intensities. Only these data enter 
the reconstruction procedure (|33fl - l|34|) . 

Let us consider the optical set-up that consist of free 
evolution to the distance d[ n , thin lens with the focal 
length /, and free evolution to the distance d ut- The lens 
has the finite aperture diameter 2a. The whole device 
can be transversely shifted to a distance s from axial 
position. The longitudinal distance d[ n from the object 
and transverse shift s are parameters of the optical set-up 
and they may be adjusted during measurement, while the 
other parameters are kept constant. The point-spread 
function h(x,xo) = h(x,xo;di u ,s) = (x\T(d- ul , s)\xq) of 
the device under consideration reads 

h(x,x ) = (x\ exp(-i^p 2 ) J\\S)®x 

k d- 
x exp (— i— x 2 ) exp Ht^P 2 ) exp (-isp) \x ), (42) 

where k = 2tt/X is the longitudinal wave number. With 
the help of the relations 1 = / da; 1 = J dp \p)(p\, 

and (x\p) = (27r) -1 / 2 exp (ixp) the point-spread function 
(I42|l can be rewritten to the form 



h(x,x ;d b 



h oc (x 1 x )£(x,x ). 



Here 



hoo(x,xo) — const e 



( a: 2 i O»0+a) 



(43) 



(44) 



is the response function of the ideal refractive focusing 
device and 



The parameter A characterizes the defocusing from the 
imaging configuration, 



dm 



1 

d ut 



1 

7' 



(46) 



the parameter 9 is related to the transverse wave number, 



x 



dn 



i-i (47) 



and the function erf() denotes the common error func- 
tion, 



erf (z) 



die" 



(48) 



In the present simulation the parameters have been 
chosen as A = 600 nm, d out = 1.5 m, / = 0.5 m, and 
a = 0.6 mm. The image (A = 0) appears at the distance 
di n = 0.75 m. However, it is blurred due to the small 
aperture. In fact, details closer than the diffraction limit 



R = C X- 



(49) 



are mapped to two spots with insufficient contrast ac- 
cording to Rayleigh's criterion. The factor C depends on 
aperture shape and coherence properties of the signal. It 
equals to 0.61 in the case of circular shape and totally 
incoherent light, or to 0.82 in the case of totally coherent 
light (Abbe's resolution limit). We set C = 0.5, what 
is equal or smaller than any classical resolution limit for 
imaging with partially coherent light. The testing ob- 
ject consists from four bright spots separated by dark 
spaces. The distance 0.15 mm between the edges of the 
central closest spots is 5-times smaller than resolution 
limit R l)49[l. The corresponding optical intensity I(q) 
is shown in Fig. ^ The off-diagonal peaks of mutual 



-2R 



-R 



R 



2R 



£(x,x ) 




(45) 



represents the correction to the finite aperture. This 
tends to unity for large aperture, linia^oo £{x, xq) = 1. 



Fig. 1. Optical intensity I(q) of the true object in the 
input plane. 

intensity T(q, q') of the testing object representing the 
cross-correlations between the spots are lower than di- 
agonal ones due to the partial coherence. The contrast 



V = (/„ 



0/(4 



of the object is V = 1. 
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The object plane is discretized by 100 equidistant 
points in the interval [—1.5,1.5] mm, or equivalcntly 
[—2R,2R]. The corresponding mutual intensity T(q,q') 
is given on the square mesh of 100 x 100 points (q m ,q' n ) 
in the process of data generation and subsequent state 
reconstruction. Similarly, the detection plane x in the 
interval [—4, 4] mm is sampled only by 64 pixels Xj for 
every longitudinal distance d ln j = (0.75 — 0.05 j) m, 
j = 0, . . . , 5, and for every transverse shift sj = (—1-2 + 
0.3/) mm, / = 0, ...,8. Using the relations (JSSJ, ffiSt , 
and (|43fl - H45l) the intensities piji = Iji(xi) in the pixels i 
can be evaluated for all the configurations of the optical 
set-up. For example, in the case of imaging axial config- 
uration (J = 0, I = 4) the intensity reveals the central 
peak with small side lobes, see Fig. |21 The correspond- 
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400 



600 



800 



iteration 



Fig. 3. The exponentially fast convergence of square dis- 
tance e during the extremization process. 




iterations, see Fig. The positions and relative intensi- 



Fig. 2. Simulated relative frequencies fi (points) affected 
by 20% of background noise sample the optical intensity 
I(x) in the detection plane (lines) for imaging axial ar- 
rangement, d m = 0.75 m, s = 0. 

ing under-sampled data /, spoiled by 20% of background 
noise are shown in the same figure. These simulated rela- 
tive frequencies fyi serve as an input for the reconstruc- 
tion procedure (|33|I - H34I) . To solve the extremal equation 
(|33|l we need to find the mutual intensity on the given 
mesh as the hermitian matrix of the dimension 100 x 100. 
As an initial iteration, the uniform incoherent superpo- 
sition of all pure states on the supposed space is used. It 
exhibits flat intensity profile. It is interesting to note that 
the final results seems to be independent of the choice 
of initial mutual intensity. In the course of repeated it- 
erations of the discretized equation 1|33[) the difference 
c = JJ dqdq'[T^ n+1 ^(q,q') - T^{q,q')} 2 between succes- 
sive iterations can be used as the criterion for terminating 
the extremization process. Numerical results show that 
the difference e reaches the level about 10 -6 after several 
tens of iterations and it reaches the level 10 -12 after ap- 
proximately 1000 iterations, see Fig. El The convergence 
improves slightly for smaller portion of background noise. 
Iterated intensity starts to reveal the four-peak structure 
after relatively small number of steps. The contrast V of 
the central part of the estimated optical intensity beyond 
the diffraction limit R reaches the value of 0.56 after 1000 




Fig. 4. Reconstructed optical intensity I(q) in the input 
plane. 

ties of bright spots in estimated object match very well 
the structure of the true object. This is demonstrated in 
Fig.0 

The numerical simulations clearly show that the pro- 
posed reconstruction algorithm is feasible and could 
be implemented. It provides considerable improve- 
ment comparing to non-statistical image processing tech- 
niques, and significantly, it yields the complete informa- 
tion in the form of correlation function. 

7. Conclusion 

The purpose of the presented paper is twofold. First 
the tight connection between wave optics and quantum 
mechanics has been emphasized. The operator language 
routinely used in quantum theory can simplify the manip- 
ulation and description of optical objects, like partially 
coherent wave and response function of optical device. 
The second goal of the contribution is the mathematical 
formulation of the reconstruction algorithm for partially 
coherent signal proceeded from the maximum-likelihood 
estimation of mixed quantum state22i2iZ4i££i2i. The solu- 
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-2R -R R 2R 



g 

Fig. 5. Contour lines (thin) of the reconstructed mutual 
intensity T(q, q'). The positions of the diagonal bright 
spots as well as the positions of off-diagonal correlations 
match the true ones (thick lines). 



tion of extremal equation by means of repeated iterations 
has been suggested. The first iteration has been explic- 
itly formulated. The proposed method never yields the 
non-physical results. 

The feasibility of the method has been verified by ex- 
tensive numerical simulations. The realistic experimen- 
tal data will be considered in the forthcoming publica- 
tion. The particular numerical example shows the good 
agreement between the true and estimated states of par- 
tially coherent light beyond the diffraction limit, despite 
of under-sampled data and 20% of background noise. 

The method is able to estimate the generic signal with- 
out any prior assumptions, utilizing only real noisy data. 
The potential applications cover wide range of optical in- 
verse problems. The method may be used for the state 
estimation of the localized mode in photonic band-gap 
structures (photonic crystals), for the determination of 
the near-field short-range correlation of the signal trans- 
mitted through random medial, and for the reconstruc- 
tion of spatial and coherence properties of light confined 
and emitted by modern laser-diode sources. Moreover, 
the general quantum origin of the method allows us to es- 
timate arbitrary continuous (discretized) partially coher- 
ent physical object described by the correlation function 
or the Wigner function. The reconstruction of de Broglie 
wave function of a particle and the optical homodyne de- 
tection of a quantum state of the light mode are typical 
examples^. In short, the method is applicable to all in- 
verse problems where the precise knowledge of partially 
coherent signal (mixed state) is essential, providing that 
the measurement device and real data are known. 
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